11  Modelos de nichos ecológicos

11.1 Introducción

Los modelos de nichos ecológicos permiten predecir la distribución potencial de especies en función de datos de presencia y variables ambientales. Puede consultar estas diapositivas sobre el tema.

11.2 Instalación y carga de paquetes

Se utilizan varios paquetes para manejar datos geoespaciales y para ejecutar la modelización.

# Paquete para acceder datos en GBIF
install.packages("rgbif")

# Paquete para acceder datos geoespaciales
install.packages("geodata")

# Paquete para mapas interactivos
install.packages("leaflet")

# Paquete para modelado de distribución de especies
install.packages("dismo")

# Colección de paquetes de Tidyverse
library(tidyverse)

# Estilos para ggplot2
library(ggthemes)

# Paletas de colores de RColorBrewer
library(RColorBrewer)

# Paletas de colores de viridis
library(viridisLite)

# Gráficos interactivos
library(plotly)

# Manejo de datos vectoriales
library(sf)

# Manejo de datos raster
library(terra)

# Manejo de datos raster
library(raster)

# Mapas interactivos
library(leaflet)

# Acceso a datos en GBIF
library(rgbif)

# Datos geoespaciales
library(geodata)

# Modelado de distribución de especies
library(dismo)

11.3 Obtención de datos de presencia

En el siguiente bloque se especifica el nombre de la especie que se va a modelar. Como ejemplo inicial se desarrolla un modelo de distribución de Bradypus variegatus (perezoso de tres dedos).

# Nombre de la especie
especie <- "Bradypus variegatus"

Para obtener los datos de presencia se emplea el paquete rgbif, el cual proporciona acceso a los datos agrupados por el Sistema Mundial de Información en Biodiversidad (GBIF).

# Consulta a GBIF
respuesta <- occ_search(
  scientificName = especie, 
  hasCoordinate = TRUE,
  hasGeospatialIssue = FALSE,
  limit = 10000
)

# Extraer datos de presencia
presencia <- respuesta$data

Opcionalmente, se puede guardar el dataframe con los datos de presencia en un archivo CSV, para no tener que repetir la consulta al API de GBIF.

# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')

El siguiente bloque debe ejecutarse si se desea leer los datos desde el archivo CSV.

# Leer los datos de presencia de un archivo CSV
presencia <- read_csv('presencia.csv')

Para manejarlo como un conjunto de datos geoespacial, el dataframe retornado por occ_search() se convierte a un objeto sf.

presencia <- st_as_sf(
  presencia,
  coords = c("decimalLongitude", "decimalLatitude"),
  remove = FALSE, # conservar las columnas de las coordenadas
  crs = 4326
)

11.3.1 Registros de presencia por país

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  ggplot(aes(x = fct_infreq(countryCode))) +
  geom_bar(
    aes(
      text = paste0(
        "Cantidad de registros de presencia: ", after_stat(count)
      )
    )    
  ) +
  ggtitle("Cantidad de registros de presencia por país") +
  xlab("País") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

11.3.2 Registros de presencia por año

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  group_by(year) |>
  summarize(n = n()) |>
  ggplot(aes(x = year, y = n)) +
  geom_line() +
  geom_point(
    aes(
      text = paste0(
        "Año: ", year, "\n",
        "Cantidad de registros: ", n
      )
    )
  ) +
  ggtitle("Cantidad de registros de presencia por año") +
  xlab("Año") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

11.3.3 Mapa de registros de presencia

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Registros de Bradypus variegatus"))

11.4 Obtención de datos climáticos

Las variables climáticas se obtienen mediante el paquete geodata, el cual brinda acceso a varios conjuntos de datos geoespaciales. En este caso, se descargan las variables bioclimáticas de WorldClim.

# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())

# Nombres de las variables climáticas
names(clima)

 [1] "wc2.1_10m_bio_1"  "wc2.1_10m_bio_2"  "wc2.1_10m_bio_3"  "wc2.1_10m_bio_4" 
 [5] "wc2.1_10m_bio_5"  "wc2.1_10m_bio_6"  "wc2.1_10m_bio_7"  "wc2.1_10m_bio_8" 
 [9] "wc2.1_10m_bio_9"  "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
[13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
[17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"

Seguidamente, se “recortan” las capas raster climáticas provenientes de WorldClim para así cubrir solamente el área en la que se encuentra presente la especie.

# Definir la extensión del área de estudio
area_estudio <- ext(
  min(presencia$decimalLongitude) - 5, 
  max(presencia$decimalLongitude) + 5,
  min(presencia$decimalLatitude) - 5, 
  max(presencia$decimalLatitude) + 5
)

# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)

11.4.1 Mapa de variables climáticas

# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
  # palette = "inferno",
  # palette = "magma",
  palette = rev(brewer.pal(11, "RdYlBu")),
  values(clima$wc2.1_10m_bio_1),
  na.color = "transparent"
)

# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
  # palette = "viridis",
  # palette = "YlGnBu",  
  palette = "Blues",
  values(clima$wc2.1_10m_bio_12),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Temperatura", "Precipitación", "Registros de Bradypus variegatus")
  ) |>
  hideGroup("Precipitación")

11.5 Modelización

Para elaborar el modelo de distribución se utiliza el paquete dismo, el cual implementa varios algoritmos de modelado de distribución de especies. Para este ejercicio, se utiliza el algoritmo Maxent. El algoritmo recibe como entradas un conjunto de datos de presencia de una especie y un conjunto de variables climáticas (en formato raster).

Primero, se eliminan las coordenadas duplicadas del conjunto de datos de presencia.

coordenadas_presencia <- data.frame(
  decimalLongitude = presencia$decimalLongitude,
  decimalLatitude = presencia$decimalLatitude
)

coordenadas_presencia <- unique(coordenadas_presencia)

Seguidamente se dividen los datos de presencia en dos subconjuntos:

• Entrenamiento: para desarrollar el modelo.
• Evaluación: para evaluar el modelo.

# Dividir los datos en entrenamiento y evaluación
set.seed(123)
n_presencia <- nrow(coordenadas_presencia)

indices_entrenamiento <- sample(
  1:n_presencia, 
  size = round(0.7 * n_presencia)
)

entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]

Ahora puede ejecutarse el modelo con Maxent.

# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima_raster <- raster::stack(clima)

# Ejecutar el modelo
modelo_maxent <- maxent(x = clima_raster, p = entrenamiento)

# Crear un raster basado en el modelo
prediccion <- predict(modelo_maxent, clima_raster)

Se evalúa el modelo.

# Extraer valores del raster de predicción para los puntos de evaluación de presencias
eval_pres <- terra::extract(
  prediccion, 
  evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)

# Generar puntos aleatorios de ausencia
ausencias <- randomPoints(mask = clima_raster, n = 1000)

# Extraer valores del raster de predicción para los puntos de ausencias
eval_aus <- terra::extract(
  prediccion, 
  ausencias
)

# Generar estadísticas de evaluación del modelo
evaluation <- evaluate(p = eval_pres, a = eval_aus)

# Graficación de la curva ROC
plot(evaluation, 'ROC')

El mapa de distribución potencial puede apreciarse en el siguiente mapa.

# Paleta de colores del modelo
colores_modelo <- colorNumeric(
  palette = c("white", "black"),
  values(prediccion),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addRasterImage( # capa raster del modelo de distribución
    prediccion,
    colors = colores_modelo,
    opacity = 0.6,
    group = "Modelo de distribución",
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>
  addLegend(
    title = "Modelo de distribución",
    values = values(prediccion),
    pal = colores_modelo,
    position = "bottomright",
    group = "Modelo de distribución"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c(
      "Temperatura",
      "Precipitación",
      "Modelo de distribución",
      "Registros de Bradypus variegatus"
    )
  ) |>
  hideGroup("Temperatura") |>
  hideGroup("Precipitación")

11.6 Recursos de interés

Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190(3-4), 231-259. https://doi.org/10.1016/j.ecolmodel.2005.03.026